A model for flexural phonon dispersion in graphite and graphene. 
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A simple model for flexural phonons in graphite (and graphene, corresponding to the limiting 
case of infinite distance between carbon planes) is proposed, in which the local dipolar moment is 
assumed to be proportional to the curvature of the carbon sheets. Explicit expressions for dispersion 
curves with full account for the long-range dipolar interaction forces are obtained and fitted to the 
experimental data using a single adjustable parameter of the model. The parameter is expected to 
depend on the ground state configuration of molecular 7r-orbitals, the same both for graphite and 
for graphene. At decreasing carbon sheet separation (high pressures) the phonon spectrum displays 
instability, corresponding to the graphite to diamond transition. Being explicitly based on the 
local dipolar moments, the proposed simple model may prove useful for considering electron-phonon 
interaction. 

PACS numbers: 81.05.U-, 63.20.-e, 63.20.dh 
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Recent interest to graphene, sparkled by its controlled 
production^ and promising electronic properties^ pro- 
duced demand for detailed study of all the related prop- 
erties of layered carbon. While conduction of graphene 
is understood relatively well, there is still a need for a 
simple model for its mechanical properties!^ In this pa- 
per such a model for out-of-plane (flexural) oscillations of 
carbon sheets is proposed with the emphasis on the long- 
range interactions between the induced electric dipoles. 
Due to orientation of the dipoles, the effect of such inter- 
actions on the flexural modes is the strongest. 

The parameters of the model have clear physical mean- 
ing and the interaction of the induced dipoles fully ac- 
counts for their mutual orientation and distance. This 
is why the model can be expected to be transferable to 
consideration of out-of-plane motion of atoms in many 
(single-, many- and few-) layered carbon allotropes, dif- 
fering only in the arrangement of atoms and, conse- 
quently, of the dipoles. As a test, the flexural phonon 
spectrum of graphite (and graphene, which is just a limit- 
ing case) is evaluated analytically and fitted to the exper- 
imental data from the literature. Despite its simplicity, 
the model correctly demonstrates an expected instability 
of graphite under pressure. 

It is well known from Chemistry that carbon has va- 
lence of four and, in its layered form, makes three strong 
chemical covalent bonds to neighbouring atoms (carbon 
or others), called the <r-bonds. The remaining electron 
also participates in bonding, forming a weak 7r-bond, fa- 
mous for making the electron delocalized across the whole 
molecule/crystal, which leads to many interesting prop- 
erties of aromatic hydrocarbons, conduction of graphite 
and metallic-like conduction of graphene. The crucial im- 
portance of 7r-bonds in hydrocarbons was realized even 
before the Hiiskel models These same bonds, as it will 
be seen from the next, define mechanical properties of 
layered carbon allotropes as well. 

The electron density around a single sp 2 -hybridised 
carbon atom is sketched as a "ball and stick" model on 
the inset in Fig. [1] . Balls represent centers of negative 




FIG. 1: Sublattices in the graphite lattice. Inset shows 
schematically the electron charge density around a single sp 2 
hybridized carbon. 



charge with er-electrons shown as —1 charged balls in 
the horizontal plane, and the 7r-electron shown as two 
— 1/2 balls above and below the plane. The thicker line 
between the —1/2 balls symbolizes electric connection 
between these charges, arising from the fact that they 
represent the same single electron. The whole picture is 
the result of momentum quantization (fixing the shape 
of electronic clouds) and simple electrostatic repulsion. 

When three more carbon atoms are connected to the 
original one, additionally to forming a bonds, their 7r 
clouds overlap, forming two "seas" of delocalized elec- 
trons above and below the plane, containing carbon nu- 
clei. 7T electrons spend half of their time above and half 
below the atom plane, and are free to move from one 
atom to the other. Evoking a-ir separability and for- 
getting about a bonds we can imagine a single layer of 
carbon as three layers of charge: a layer of +1 (per atom) 
charges, representing the uncompensated charge of car- 
bon ions, and two -1/2 charged layers of ir electrons on 
both sides of it. 

Having this picture in mind, imagine that such a tri- 
layer carbon sheet is curved. Because both "seas" of elec- 
trons are connected (it is the very same electrons after 
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all) the charge, pushed by the electrostatic repulsion, is 
free to redistribute from the contracting to the expanding 
side. This creates a local electric dipolar moment (inter- 
acting with similar dipolar moments across the layer) , in- 
creases electrostatic energy of deformed electronic clouds 
and of the layer as a whole, producing the restoring force. 

To describe this process mathematically, consider a 
graphite-like lattice, split into four independent hexag- 
onal Bravais sublattices 1-4, shown in Fig. [TJ These 
sublattices are essentially the same, but only shifted 
with respect to each other, so that position of a lat- 
tice site, identified by three-dimensional integer vector 
rh = k] € Z 3 and number of the sublattice I, is 
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where a is the nearest neighbour distance in the lattice 
planes, a is dimensionless inter-plane distance (in units of 
a), u~ is out-of-plane displacement of atoms (in units of 
a), the matrix (denoted in the further text as b) contains 
basis vectors of the lattice. The sublattice displacements 
(in units of a) are 



d x = 



(2) 

All neighbours of an atom at sublattice I belong to an- 
other sublattice I (by definition: 1 = 2, 2 = 1, 3 = 4, 
4 = 3). For the lattice ([T]) three nearest neighbours 
of an atom rh on sublattice I are the atoms rh, rn^ = 
rh + A;[l, 0, 0], and rhf = rh + A;[0, 1, 0] of sublattice I, 
where A; = sign(l — I); sign(x) is 1 if x > 0, —1 if x < 0. 

There can be several definitions of local surface cur- 
vature, but, for the case of small deformations of the 
original lattice, all of them are essentially the same up 
to a constant multiplier. It is convenient to measure the 
curvature as a distance of the considered atom from the 
plane, defined by its three nearest neighbours. The nor- 
mal to this plane at site rh of sublattice I is proportional 
to 
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where cross denotes the vector product. The local dipolar 
moment is then proportional to 



I rh 
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where dot stands for the scalar product. Up to the first 
order in atom displacements u we get p — [0, 0,p], where 
p is 
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The Hamiltonian is then 
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where m is an atom's mass, c and b are parameters of 
the model (both have dimensions of energy). Express- 
ing this Hamiltonian in units of ma 2 , we can introduce 
two characteristic frequencies ujq=c/ (ma 2 ) and u± — /3u>q 
with ft = b/c. The parameter ujq defines the overall en- 
ergy scale (later we normalize it out), while j3 remains 
the free parameter of the model. 

Physically, the model attempts to capture essentials 
of 7r-orbitals polarization during the deformation of each 
individual graphene sheet. Such deformation produces 
local shift of the charge from one side of the sheet to an- 
other, which can be represented as an additional charge 
density, superimposed over the original, undeformed, or- 
bital. The first potential energy term in Q corresponds 
to the electrostatic self-energy of this additional density, 
while the second term models the interaction between 
these redistributed charges across the whole lattice. This 
reproduces precisely the extremely short-range (the self- 
energy, taken simply as an independent parameter) and 
long-range parts (by keeping the leading-order dipolar 
terms) of the interaction between deformed orbitals while 
neglecting the higher-order multipole terms, whose con- 
tribution peaks at intermediate distances. 

To solve the model one may reexpress the Hamiltonian 
© in terms of the displacements uL and differentiate to 
find the force on an element rh of each of the four sub- 
lattices. Representing the displacements by their Fourier 
components both in time and space 



i4(t)= / U '(fc)e 2 ™^+^d 3 £, 
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where k = [kx,kY,kz] and the explicit dependence on 
time t is shown, one gets the following usual secular equa- 
tion for the frequency 
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where the matrix (called the dynamical matrix and de- 
noted here, including the numerical coefficient 1/9 in 
front, as M) is obviously self-adjoint. Its elements are 
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= 2a (2 + pSi) - 3^ (bS 2 + bS 2 ) 
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FIG. 2: Dispersion of flexural phonons in graphite across the 
first Brillouin zone, labels correspond to the well known high- 
symmetry points. Experimental data are shown by circles, 8 
squares^ (as reproduced in Reffl3). triangles, 10 diamonds.— 
Solid lines are calculated from (|8))- (|14[) with ft = 0.36, dashed 
lines are corrected by (|15p with r\ = 0.036. 



with Si = Z b (3,k,di), Ui = Si - 3cv 2 Z b (5, k, di), a = 
2 cos (V3nkx) cos(37rfcy) + cos (2^/3nkx) + 6 and b 
I + 2e~ mkY cos (%/37rfc x )- This assumes the following 



definition of the Epstein zeta function 
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where prime near the sum means that singular terms are 
excluded, A is an arbitrary DxD matrix, s is (in general) 
complex scalar and c, d are arbitrary Z?-vectors. The 
vectors di are from ([2]) with do = [0, 0, a]. 

Epstein zeta function can be very efficiently evaluated- 
by a computer program^ The four branches of flexural 
phonon spectrum of graphite at any point in /c-space are 
then just the square roots of eigenvalues of the matrix M, 
defined by (|5|)- p^|) . Apart from the parameter u>o, defin- 
ing the overall frequency scale, these branches depend on 
graphite interlayer separation a = 2.34, taken from the 
experiment, and the free parameter of the model ft. 

A wealth of experimental data on phonon spectrum of 
graphite is available in the literature^— Some of these 
data are shown in Fig. [2] along with dispersion curves 
predicted by the considered model for ft = 0.36, giving 
the best fit to the data. The value of loq = 51meV was 
obtained by fixing the values of the spectrum at K point. 
There is a slight disagreement near M, which can be at- 
tributed either to the oversimplification of the model, ne- 
glecting the higher-order multipolc terms, or may be even 
to the experimental errors, which, at least for one set of 
measurements^ increase on approach to M. Provided 
there is a single adjustable parameter, the agreement be- 
tween the model (solid line) and the experiment is very 
good. 

The spectrum of graphene can be obtained as a limit 



at a — > oo. Then C,D,E (interaction between 
the layers vanishes) and M splits into two 2x2 sub- 
matrices. Expressions for A and B remain the same, 
except that zeta function becomes 2-dimensional as the 
matrix b in the expression for Si loses its last row and 
last column. The resulting spectrum is very similar to 
the one already shown in Fig. [21 just there is no split- 
ting of acoustic and optical branches. Please note that 
even though the long-range interaction across the layers 
is eliminated in graphene limit, the interaction inside the 
layer still remains, giving the spectrum its specific shape. 

Let us also note that in the present model displacement 
of graphite layers as a whole, without flexing them, does 
not generate any dipolar moment ^ and, thus, leaves 
the energy invariant. According to the Goldstone the- 
orem, the presence of such continuos symmetry implies 
the existence of an additional acoustic mode in the spec- 
trum (the other acoustic mode is due to the energy in- 
variance with respect to translation of the whole crystal). 
In graphite this mode turns optic, acquiring a certain 
amount of energy at T point due to macroscopic van der 
Waals interaction between the carbon layers and depends 
on their complete flexural phonon spectrum (|5])- (fl~4")) as 
well as temperature. This interaction can be simulated 
by introducing an additional phenomenological harmonic 
coupling between neighbouring atoms on sublattices 2 
and 3, resulting in the following addition to the dynamic 
matrix in (El) 
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where rj is a parameter. The effect of this addition is 
mostly localized in the neighbourhood of T point (see 
dashed line in Fig. [51 corresponding to the best-fit value 
of r) = 0.036). In principle, the value of rj can be calcu- 
lated on the basis of the present model (by also including 
the repulsion between the layers due to exchange interac- 
tion), but, while its introduction results in a better fit to 
the data, it is methodologically wrong to try to model the 
macroscopic interaction inside the microscopic Hamilto- 
nian 

One may try to add er-bonds stretching to the model 
by introducing the following term into the Hamiltonian 
©: 
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where 7 is a free parameter. This models a-bonds as 
harmonic strings with an equilibrium length of 1, produc- 
ing, up to the first order in displacements, an additional 
force —672?^ at each site (proportionality of this force 
to the magnitude of the local dipolar moment is just a 
convenient coincidence). The corresponding dynamical 
matrix can be obtained by adding 37 to A in Eq. @ 
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and -7(1 + e z " kY cosy/Znkx) to B in Eq. JTDJ. How- 
ever, fitting the resulting spectrum to the experimental 
data, produces (up to the fitting error) 7 = 0. It means, 
there is no contribution of er-bonds stretching to flexural 
phonon spectrum, which may be a consqeuence of their 
strong anisotropy. There is no need to introduce an ad- 
ditional parameter 7 to obtain the fit, shown in Fig. [__ 

The interlayer interaction becomes progressively 
stronger if one presses the graphite, reducing the inter- 
layer distance a. In this case splitting of the acoustic 
and optical branches rapidly grows, until, at a critical 
value of a = 0.91, one of the acoustic branches touches 
the horizontal axis at point K. At smaller a the cor- 
responding eigenvalue of M becomes negative, meaning 
that the lattice is unstable. This happens at interlayer 
distance approaching the intra-layer distance between 
carbon atoms, suggesing that it corresponds to lability 
boundary of graphite — > diamond transition. 

To conclude, the presented simple model, by explicitly 
including the long-range dipolar interactions, quantita- 
tively reproduces flexural phonon spectrum of graphite 
and graphene using a single, universal for all layered 
carbon allotropes, parameter /?. This is an advantage 



with respect to the widely used for this task Born-von 
Karman type models, necessitating to include several 
nearest neighbours (and, consequently, many parame- 
ters) into consideration to obtain comparable agreement 
to the experiment. Because the interactions in the pre- 
sented model are purely electrostatic, their dependence 
on inter-atomic and inter-layer distance is explicit. The 
model is not specific to graphene or graphite. The ex- 
pression for stress-induced dipolar moment and a similar 
Hamiltonian may be useful while considering other lay- 
ered carbon allotropes, such as single- and multi- walled 
nano-tubes, fullerenes etc. 

The effect of higher order multipoles or other quickly 
decaying short-range interactions would correspond to a 
certain simple addition to the dynamical matrix, simi- 
lar to discussed above. Since the most important (and 
more difficult to evaluate) long-range contribution is al- 
ready taken into account (and produces a very complete- 
looking spectrum), one can expect the other contribu- 
tions to be small and local. Thus, it can be hoped that 
the model can be a solid basis for further quantitative 
improvement. 
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